Linking physics and algorithms in the random-field Ising model 
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The energy landscape for the random-field Ising model (RFIM) is complex, yet algorithms such as the push- 
relabel algorithm exist for computing the exact ground state of an RFIM sample in time polynomial in the 
sample volume. Simulations were carried out to investigate the scaling properties of the push-relabel algorithm. 
The time evolution of the algorithm was studied along with the statistics of an auxiliary potential field. At very 
small random fields, the algorithm dynamics are closely related to the dynamics of two-species annihilation, 
consistent with fractal statistics for the distribution of minima in the potential ("height"). For d = 1,2, a 
correlation length diverging at zero disorder sets a cutoff scale for the magnitude of the height field; our results 
are most consistent with a power-law correction to the exponential scaling of the correlation length with disorder 
in d = 2. Near the ferromagnetic-paramagnetic transition in d = 3, the time to find a solution diverges with a 
dynamic critical exponent of z — 0.93 ± 0.06 for a priority queue version and z — 0.43 ± 0.06 for a first-in 
first-out queue version of the algorithm. The links between the evolution of auxiliary fields in algorithmic time 
and the static physical properties of the RFIM ground state provide insight into the physics of the RFIM and a 
better understanding of how the algorithm functions. 



I. INTRODUCTION 

Models of materials with quenched disorder, such as spin 
glasses and random magnets [1], typically have extremely 
slow dynamics in the low-temperature glassy phase, due to the 
existence of many metastable states separated by barriers that 
grow with the system size. Such exponential slowing down af- 
fects optimization methods that are modeled on the dynamics 
of the physical system, such as simulated annealing Q]. The 
slow dynamics of the system prevent precise study of the equi- 
librium behavior. Besides using physically motivated dynam- 
ics, one can consider using combinatorial methods 1 3||4j] to 
compute partition functions or ground states. Some problems, 
such as determining the ground state of a 3D spin glass |5] or 
finding the partition function for the RFIM at finite temper- 
ature, are NP-hard |6, 7]. Though combinatorial approaches 
derived in computer science greatly accelerate searches for the 
ground state (HEl], it is difficult to study systems with more 
than 10 3 degrees of freedom. However, many problems for 
disordered materials, such as computing the partition function 
of a 2D spin glass |5] or finding the ground state of an RFIM 
sample 1 10], can be solved in time polynomial in the volume 
of the sample, allowing for the solution of samples with over 
10 7 degrees of freedom. 

The push-relabel (PR) algorithm introduced by Goldberg 
and Tarjan |11] is directly applicable to finding the ex- 
act ground state of the RFIM and has been extensively ap- 
plied to study the RFIM's zero-temperat ure p aramagnetic- 
ferromagnetic transition Q2 [lj H [H Qjf . The 
running time is bound by a polynomial in the number of spins, 
but there is a power-law critical slowing down of the PR algo- 
rithm at the zero-temperature (T = 0) transition fl3lflrifl9il . 

This paper explores the connection between the auxil- 
iary variables of the PR algorithm and the zero-temperature 
disorder-driven phase transition. We also look at the algorithm 
in the limit of small disorder, where the dynamics of the algo- 
rithm turn out to be closely related to annihilation processes 
studied in statistical physics |20, 

num. 

We define the RFIM Hamiltonian, review its phases, and 



define the rules for the push-relabel dynamics in Sec. [[J These 
rules describe the evolution of auxiliary fields; the dynamics 
of these fields leads directly to a ground state for the RFIM. 
Our main focus will be on the "height" or potential field. This 
field guides the determination of the spin-up and spin-down 
domains in the ground state. Given these definitions, we out- 
line our results in Sec. [II]] We study the relationship of the PR 
algorithm to annihilation processes at low values of the disor- 
der by studying the algorithm dynamics and final topography 
of this auxiliary potential surface in Sec. II VI In Sec. [V] we 
study the topography of the potential surface for general dis- 
order, especially near the paramagnetic-ferromagnetic transi- 
tion. The primary results of the paper, especially the scaling 
of the running time and statistics of the potential field, are 
summarized in Sec. lVIl 



II. RANDOM-FIELD ISING MODEL AND THE 
PUSH-RELABEL ALGORITHM 

The random-field Ising model (see, e.g., 1 1] and references 
therein) captures essential features of models in statistical 
physics that are controlled by disorder and have "frustration", 
i.e., energy competition between different terms of the Hamil- 
tonian. Such systems have complex energy landscapes and 
the ground-state configuration for a given sample is not usu- 
ally obvious. Large barriers separate a very large number of 
metastable states. If such models are studied using simula- 
tions mimicking the local dynamics of physical processes, it 
takes an extremely long time (exponential in a power of the 
system size) to encounter the exact ground state. But in many 
cases there are very efficient methods for finding the ground 
state. These methods break away from a direct physical rep- 
resentation. Extra degrees of freedom are introduced and an 
expanded problem is solved. By expanding the configura- 
tion space and choosing proper dynamics, the algorithm goes 
"downhill" in a fashion that avoids having to go over barriers 
that exist in the original physical configuration space. An at- 
tractor state in the extended space is found in time polynomial 
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in the size of the system. When the algorithm is completed 
by finding this attractor or minimum in the extended space, 
the auxiliary fields can be projected onto a physical configu- 
ration, which is the guaranteed ground state. The RFIM is an 
example where this extension can be carried out. 

A. Model 

In the RFIM, there is competition between ferromagnetic 
terms, characterized by a strength J and local random fields 
hi of characteristic magnitude J A, with a Hamiltonian 

h = - jy] sjSj - hjSj, (i) 

(y) 1 

where the Ising spins Sj = ±1 lie on sites i on a d-dimensional 
lattice and the notation (ij) indicates a sum over nearest- 
neighbor pairs of sites. The samples have linear dimension 
L, with n = L d sites, and we use periodic boundary condi- 
tions. The independent Gaussian random variables hi have 
mean and variance A 2 J 2 . 

In dimensions d > 2, there is a zero-temperature transi- 
tion between two phases at the critical disorder A = A c . 
When A < A c , the ferromagnetic interaction between near- 
est neighbors dominates and the spins take on a mean value 
to = n^ 1 s i w i m \ m \ 7^ in the limit n — > oo. In the 
case A > A c , randomness dominates and the ground state is 
"paramagnetic", with \m\ = 0, as n — > oo. In the standard 
picture, the zero-temperature transition has the same critical 
exponents as the finite temperature transition. 

In dimensions d = 1,2, there is no zero-temperature tran- 
sition. For a given sample size L, there is a characteris- 
tic crossover value A X (L) for the disorder strength. When 
A <C A X (L), samples of size L have a high probability to 
have uniform spin. For large disorder, A 3> A X (L), there are 
many domains of uniform spin in the ground state. Exact cal- 
culations and scaling arguments show that A X (L) ~ L -1 / 2 
for d = 1, while scaling arguments and computer simulations 
give A X (L) ~ [ ln(L)1~ 1//2 for d = 2, up to logarithmic cor- 
rections ^iHlHEiElSlii. 



B. Push-Relabel Algorithm 

We present here a description of the auxiliary fields and al- 
gorithmic dynamics used in the push-relabel algorithm. We 
do not provide the standard proof of the correctness of the 
algorithm (see 0, H [3 l30ft for such proofs and further dis- 
cussion). To describe our results, it is necessary to define the 
auxiliary variables used and the update rules. It is also use- 
ful to provide an intuitive description of the dynamics of the 
algorithm. 

There are three auxiliary fields: (i) the excess e*, (ii) the 
residual strength (capacity) defined for ferromagnetically- 
coupled pairs (ij), and (iii) a distance or height field it,-. Ini- 
tially, the fields are set according to the rules = hi, r,j = J 
and Ui = 0. A site i is "active" if e, > and Ui < oo; a site is 



a "sink" if ej < 0. Primarily two types of operations, "push" 
and "relabel", modify the fields at active sites and their neigh- 
bors. The push operation rearranges the locally conserved 
and also modifies the . The conditions for a push from site i 
to neighboring site j are that U{ = Uj + 1, ej > and > 0. 
If a push is executed, the excess e, is modified via ej — > ei—Si, 
where Si — min(ej, n 3 ■), and the residual bond strengths are 
modified according to r,j — > r,j — Si, and ry L — > rj% + The 
relabel operation increases the it, at an active site i. Whenever 
a push at an active site i is not possible, m is increased to the 
minimum value to enable a push from i (Ui is set m — oo, if 
no push is possible from i due to saturated bonds). As seen 
from the rules for a push, the height field U{ guides the rear- 
rangement of excess. The pushes are always downhill with 
respect to the height field. 

As a motivation, though again, not a proof, pushing the ex- 
cess field roughly corresponds to a rearrangement of the ex- 
ternal magnetic fields. Regions that have large ferromagnetic 
coupling compared to A will have a net excess that is of the 
same sign as the total magnetic field over the region. The 
pushes allow the spin orientations in these domains to be de- 
termined by cancellation of positive and negative excess. Note 
that as pushes are allowed only when > 0, the rearrange- 
ment of excess is limited. This leads to finite size domains, 
for large enough samples and strong enough disorder. 

The set of rules given so far does not define an algorithm. 
In order to have a well-defined procedure, one must organize 
the push and relabel operations and set a criteria for termina- 
tion. When considering a site with > 0, we first execute 
all possible pushes and then relabel site i, if necessary. This 
is defined as a single push-relabel step; the number of such 
steps will be our measure of algorithmic time. The order in 
which sites are considered is given by a queue. In this paper, 
we compare results for two types of queues: a first-in-first- 
out queue (FIFO) and a lowest height priority queue (LPQ) 
(other options are available 1171 Oil 13211 ). The FIFO structure 
executes a PR step for the site i at the front of a list. If any 
neighboring site is made active by the PR step, it is added to 
the end of the list. If i is still active after the PR step, it is also 
added to the end of the list. This structure maintains and cy- 
cles through the set of active sites. The LPQ structure is also 
a list, but is always sorted by the height label it,-. The sites 
with lowest it, are always at the front of the list. The LPQ 
version executes PR steps for active sites that are likely to be 
near sinks. This algorithm will act repeatedly on the set of 
sites with lowest height. 

The algorithm terminates when no active sites remain. Sites 
with positive excess will remain, in general, but the height at 
those sites will be u, = oo. At the end of the algorithm, the 
sign of spin Si is set to be positive if there is no path from i to 
a site with negative excess along bonds (uv) with r uv > 0. If 
there is such a path, the spin Si = — 1 in the ground state. 

One operation that greatly speeds up the running time of 
the algorithm is the global update 14 fill [171 13011 . This oper- 
ation recomputes the height field u, so that it, is the minimal 
distance from the site i to the set of sinks. The height is set 
to Ui — oo when there is no path from i to a sink along with 
all edges satisfying r uv > 0. We denote the period between 
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global updates as T. In the rest of this paper, we fix T = n 
(d = 2, 3) or r = 2n (d — 1), which gives near minimal 
running times for the PR algorithm with fixed global update 
intervals OIH. 



III. OVERVIEW OF RESULTS 

Our principle results concern the statistics of the height field 
Ui and the connection between the topography and the run- 
ning time of the ground state algorithm. The local relabels 
and global updates result in a height field that guides posi- 
tive excess to maximal cancellation with the negative excess 
sinks, subject to constraints from the residual bond strengths 
Tij . The final configuration of the height field ui determines 
the ground state, though a single ground state can be consis- 
tent with a large set of terminal height field configurations. 
The choice of the order of operations will determine the final 
height field in a given sample. The time scale needed to estab- 
lish this height field, while cancelling excess fields and modi- 
fying residual bond strengths, determines the running time of 
the algorithm. In this paper, we mostly study the FIFO data 
structure. Besides being the fastest data structure we used, this 
structure is most natural for making connections to physical 
dynamics, where particles move in parallel. For this structure 
and algorithm, sites with positive excess are each treated once 
during a cycle through the active sites. This is to be contrasted 
with a data structure like LPQ, where a single parcel of excess 
may be moved many times while other portions of the system 
remain static. 

Using the FIFO data structure, we first study the limit of 
small A, where the rearrangement of excess is not affected by 
the bond strength. By rearranging the positive excess, the al- 
gorithm cancels out negative and positive excess as much as 
possible. We arrive at the natural connection that the dynam- 
ics of the algorithm at weak fields is related to the extensively 
studied set of annihilation models A + B — > 0. In such mod- 
els, there are two types of particles, A and B, with one or both 
types mobile. The particles annihilate (or combine to an inert 
particle) upon contact. In general, the motion of the particles 
is modeled as due to random diffusion or to overdamped drift 
caused by interactions between the particles. A very rich set 
of scaling results have been found to describe the dynamics 
of the average density, domain sizes, and domain profiles in 
the annihilation process lEol EH E2I l33l l34ll . For a descrip- 
tion of the push-relabel algorithm, we can consider positive 
excess as A particles and the sinks as B particles. The B 
particles are immobile. For Gaussian disorder, the particles 
will not exactly annihilate: upon meeting, either the positive 
excess or negative excess will be saturated by the excess of 
opposite sign. However, the cancellation of negative and pos- 
itive excess leads to a decrease in the density of active sites 
similar to the direct cancellation A + B — > 0. Particles of 
type A may coalesce, changing the speed at which the densi- 
ties evolve. In the small disorder limit, we find that the run- 
ning time grows very slowly - apparently logarithmically in 
L - when measured in the number of PR steps per site. The 
final distribution of sinks is found to have a fractal charac- 



ter at small length scales; this fractal character is consistent 
with that for annihilation processes |21], at least when d = 1. 
This fractal distribution is related to a power-law distribution 
for the height values. We also study the time evolution for 
the density of A and B particles. When the A particles can 
coalesce, the result is a single excess packet of large weight 
at long times. Forbidding this coalescence by modifying the 
algorithm leads to an approximate equal density of sinks and 
active packets during the solution process. 

We also studied the number 2Vpr of push-relabel operations 
required to find the ground state and the topography of the 
height field Ui for general A. For d = 1,2, the peak running 
times were used to define the crossover field A x . The size 
dependence of A^ is consistent with the expectations A^ ~ 
L- 1 ' 2 (d = 1) and L ~ (A x )~ 2y e~ A o/ A *, with y « 1 
and a fitted value for Ag w 1.3. The scaling of the height 
fields in d = 1, 2 is consistent with this same divergence in 
correlation length £: the fraction of sites P(u) with height u 
is proportional to a function of 

For d = 3, the running times near A c exhibit distinct dy- 
namic critical exponents for the running times. Excellent scal- 
ing results when we use values for the critical disorder A c 
and correlation length exponent v derived from more physical 
measures of the RFIM ground state JkJ- F° r LPQ we find 
the dynamical critical exponent z = 0.93 ± 0.06, while for 
FIFO, z = 0.43 ± 0.06, where z describes the running time 
via Npr ~ L z . The probability distribution for height fields 
decreases exponentially with height u when A > A c , while 
it is increasing at small u for A < A c . At the critical point 
A = A c , the probability distribution P(u) is very nearly con- 
stant out to the linear size of the system. We also study the 
structure of the domains by analyzing the paths to the sinks 
near A c ; these paths are apparently nonfractal for all A. 



IV. SMALL DISORDER LIMIT 

The relationship between the auxiliary fields and running 
times is simplest when the random field is weak compared 
to the magnitude of the exchange coupling. In this limit, the 
residual capacities of the directed links connecting sites are 
never saturated by pushes and, as a result, the push operations 
are unrestricted. The dynamics, then, is roughly described by 
the motion of positive excess towards the sinks. This would 
be exactly true if global updates were carried out immediately 
whenever a sink was removed by annihilation with positive 
excess. The height field would then always guide positive ex- 
cess directly towards the nearest negative sinks. Note that to 
maintain this exact guidance one need do even less: upon an- 
nihilation of a sink, the height field only in the region that 
acted as a funnel for that sink would need to be updated. We 
restrict ourselves to the standard approach using only local re- 
labels and global updates at periodic intervals. 

For much of the evolution time of the algorithm, the global 
updates maintains a m landscape that approximates well the 
one that would be found for more rapid updates (at least at 
early times, when there are many active sites). We note that 
for very weak disorder, where residual bonds are never sat- 
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urated, a global update constructs a height field that exactly 
equals a potential field where the negative-excess sites (the 
sinks) act as sources. This potential is the minimum over all 
sinks of a potential that increases linearly with distance from a 
sink. The packets of positive excess at active sites are guided 
by this potential from the sinks but do not interact with each 
other via any potential. This lack of repulsion between the 
A particles is one difference from the force-guided motion 
for annihilation processes that has been previously considered 
lEoll . As already noted, the positive excess packets do interact 
by coalescence when an excess is pushed onto a site already 
containing excess. 

Numerically, we studied the low disorder limit using two 
varieties of the algorithm. We first varied A and examined the 
A — > limit. To determine what disorder parameters gave 
this limit, we examined a wide range of values with A <C 1. 
We used this data to determine a value of A where the quan- 
tity of interest (such as Npr) was constant over a factor of at 
least 10 3 in A. This assured that the small-disorder limit had 
been reached, but that A was not so small that the discrete- 
ness of the disorder distribution affected the running time. In 
our algorithm, integer values for J and hi were used, with 
J = 5 x 10 8 . The product J A characterizes the integer reso- 
lution for the disorder; small values do affect the results. For 
example, a small decrease in Npp was seen for J A < 10 2 . 
In samples of less than 2 x 10 6 spins, we found a range of A 
where the running time and other quantities were quite con- 
stant. For the running time data reported here, where finite J 
is used, a value of J A of value 10 3 or 10 4 was well within this 
range. Some quantities, such as the number of sites with pos- 
itive excess at the end of the algorithm, required even smaller 
values of A at fixed J, in the largest samples. This may be 
somewhat surprising, as the magnitude of the excesses that 
are rearranged is still always much less than J, so that a single 
packet of positive excess will not saturate a bond. However, 
some bonds end up being on the path of many packets, so that 
their residual bond strength is driven towards zero. Because 
of this detail, we compared our results against a second pro- 
gram in which J = oo. In this version of the algorithm, there 
were no limits set on the capacity of a bond. Besides allowing 
one to study the A — » limit directly, this code is also sim- 
pler than the full push-relabel code and requires less memory, 
as we do not need to maintain the 2d values of ry at each 
site, allowing us to study samples up to size 512 3 . We ver- 
ified the results of this simpler code on a sample-by-sample 
basis in many cases to verify that the results agreed with the 
A — > limit of the full push-relabel program which did not 
take J — > oo. After confirming the correctness of the newer 
approach, we used it to generate most of the data used in this 
subsection. 



A. Running time at small A 

A direct way to measure the dynamics of the algorithm is 
to examine the dependence of the running time, measured by 
the number of push-relabel operations iVpR, on system size 
L. We first present such data for the case of very weak disor- 



der. The dependence of Npp/n = NppL~ d on L is plotted 
for d = 1, 2, 3 in Fig. ^ F° r the ID systems, we studied 
samples over the size range L = 2 — > 2 24 . The data is con- 
sistent with an asymptotic approach to L -1 JVpR ~ ln(i), 
though the apparent slope of the L _1 iVpR vs. L plot becomes 
approximately constant only for larger L > 5 x 10 3 . For the 
last two decades in scale, a logarithmic fit describes the data 
well. The growth in iVpR with L for d — 2 is also very slow 
compared to all but the smallest power laws. The slope on a 
linear-log plot of Npji/n vs. L is not constant over the length 
scales studied, but the data is not inconsistent with conver- 
gence to Npr ~ I? ln(L) for L larger than a crossover point 
L x ~ 10 3 , similar to the d = 1 crossover in shape and at 
roughly the same Np^/n or L. The 3D data is also consis- 
tent with very slow growth for the running time at small A, 
but given the d = 1 and d = 2 results and that we are re- 
stricted to smaller L < 512 and Npp/n < 5, it is not pos- 
sible to make a strong conclusion on a functional form for 
Npp,(L, A — > 0) when d = 3. Over the size range that we 
studied in the case d = 3, the slope is not smoothly varying 
on a linear-logarithmic plot of Np^/n vs. L, due to the dis- 
creteness resulting from global updates seen at small Np^/n, 
but the approximate behavior is logarithmic over this range of 
scales. 

A plausibility argument can be made for a logarithmic de- 
pendence of iVpft/n on L. If the algorithm successively 
"solves" for the ground state by coalescing positive and nega- 
tive excesses, unhindered by residual capacity constraints, one 
might expect that the coalescence is carried out on succes- 
sively larger scales, leading to a ln(X) number of intermediate 
solutions. The pushes lead to cancellation of excess on a given 
scale. If each length scale requires a constant number of push- 
relabel operations per site, this would give Npr ~ L d ln(L). 
In d = 1, this is the most likely scenario, as rearranging ex- 
cess over a scale £ requires a minimum of £ PR steps. In higher 
dimensions, however, not all sites must be activated in a vol- 
ume £ d for the positive excess to rearrange on a scale I. A 
better understanding of the dynamics than is presented here is 
needed to confirm this tentative description in higher dimen- 
sions. To explore the dynamics further, we next consider how 
many sites with positive excess are being rearranged and the 
number of sinks present. 



B. Number of remnant sinks and sources 

The ground-state magnetization is uniform but has random 
sign in a finite system in the limit of weak disorder. The topog- 
raphy of the height field, however, is very different between 
the two possible ground states (up and down). If the spins are 
all positive at the completion of the push-relabel algorithm, 
there are no sites with negative excess, and itj = oo at all i, as 
all sites are "unreachable". If the magnetization is uniformly 
negative, there are no sites with positive excess, Ui is finite at 
all i, and there is a remnant set of negative excess sites, i.e., 
the sinks that have survived annihilation up to the completion 
of the PR algorithm. The number of sites with positive and 
negative excess gives some indication of the dynamics of the 
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Figure 1 : Dependence of the running time, measured by the number 
of push-relabel steps per spin, NpB./n, on the system size, in the 
limit of small disorder A, for dimensions (a) d = 1 (L = 2 — * 2 24 ), 
(b) d = 2 (L = 2 -> 4096), and (c) d = 3 (L = 2 -> 512). The 
FIFO data structure was used in each case. The error bars indicate 
the la error in the statistic; the minimum number of samples at the 
largest sizes is 1200, 3700, and 876 for d — 1, 2, and 3, respectively, 
with more samples at smaller L. The dependence is suggestive of a 
logarithmic dependence JVpr. ~ L d \n(L), especially when d = 1, 
as indicated by the fit in (a) for 2 24 > L > 2 18 . 



Figure 2: [Color online] Plot of L~ d/2 N-(L), the average num- 
ber of sites with negative excess at the termination of the algo- 
rithm scaled by the expected average L d//2 , vs. system size L, for 
d — 1,2,3 samples with negative magnetization, in the A — » 
limit. The plot indicates convergence to a single value for the scaled 
variable in each dimension, consistent with iV~ ~ L d ' 2 . 



algorithm. 

As the hi are independent variables with identical distri- 
bution, the magnitude of the sum of the hi scales simply as 
~ L d / 2 JA. In the FIFO and LPQ algorithms, the negative ex- 
cess is not mobile and so cannot coalesce. If the mean negative 
excess at a remnant sink at the end of the algorithm is of order 
- J A, as it is at the beginning of the algorithm, we expect that 
~ L d l 2 sites with negative excess remain. This is the most 
natural assumption, taking the cancellation of negative excess 
to be with packets of positive excess that are either compa- 
rable or much larger in magnitude. To verify this assump- 
tion, we may simply count the number of sites with negative 
excess at the end of the algorithm in samples with negative 
magnetization. We will refer to the average of this quantity as 
N~(A, L). Plots of L- d ' 2 N~ ( A = Q,L) = L- d / 2 N-{L) 
for d = 1,2.3 are displayed in Fig. [5] In all dimensions stud- 
ied, there is a convergence to a single value. This is consistent 
with the picture that the mean negative excess at a sink con- 
verges to a single value ~ — J A as L — » oo. 

In contrast, for the samples with positive magnetization in 
the ground state, the number of sites N + (L) with positive ex- 
cess does not vary rapidly with L. The number of packets of 
positive excess at the termination of the algorithm is generally 
small, with the average of N + (L) < 10 in all sizes L and 
dimension d that we examined (see Fig.|3jl. As the net posi- 
tive excess is the sum of hi in these samples, the amount of 
excess per site is large, w 0{L d / 2 A). The initial excesses of 
magnitude A coalesce significantly under the push operations 
at small A. At initial times and high density of sites with pos- 
itive excess, this coalescence may take place due to "chance" 
collisions that depend on the order of the operations at the 
sites. At lower densities, the coalescence must result from a 
focusing caused by the topography of the height landscape u.- L . 

We also studied how densities of sinks and active sites con- 
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Figure 3: [Color online] Plot of N + (L), the number of remnant sites 
with positive excess in samples with positive magnetization, vs. L 
for d = 1,2, 3, in the limit A — > 0. The small numbers indicate 
that most of the positive excess is collected into less than 10 sites 
by the termination of the algorithm. The ID data appear to converge 



to a constant N + {oo) = 5.96 ± 0.03, 
apparently converge to N + (oo) = 1 



while the 2D and 3D data 



verge to their final values N~L~ d and N + L~ d . We have 
computed the number densities p ± (npR,L) of the positive 
and negative sites as a function of "time" npR (number of 
PR cycles completed so far; iVpR = iipr at the end of the al- 
gorithm). The data, plotted for d = 2 in Fig.|4] shows that the 
number density of sites with positive excess reaches a plateau 
consistent with constant p + = L~ d /2 in the larger samples, 
i.e., N + m 1 in samples with positive magnetization. In the 
late-time regime of few packets (that is, few active sites), the 
positive excess packets may move across the sample many 
times between global relabellings, in d = 2, 3. The density 
of negative excess sites decreases rapidly with the number of 
PR operations per site, though more slowly after the positive 
excess has coalesced. 



C. Spatial structure of the remnant sinks 

We can study the topography of the Ui in the half of the 
ground states that have negative magnetization at small A. 
One approach is to study the probability distribution and cor- 
relations of Ui at all sites. We carry this out for d = 1. Very 
closely related information can be found from the locations of 
the sites where u. L = at the termination of the algorithm. The 
spatial distribution of these remnant sinks reflects the history 
of the cancellation process. We have examined this spatial 
distribution in d = 1, 2, 3. The terminal height field Ui can 
also be computed from the final sink locations alone, when 
A — * 0, so these two descriptions are closely related. One 
method to study the distribution of the sinks is to coarse-grain 
the distribution of the remnant sinks over a length b. The com- 
puted quantity we used, N(b), is defined as the number of 
non-overlapping boxes of dimension b d that contain at least 
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Figure 4: Plots of the number densities p + and p~ vs. jipr, the num- 
ber of PR steps executed during the algorithm for (a) positive excess 
packets (active sites) and (b) sinks, for d — 2 FIFO, in the limit of 
small disorder. The density of positive sinks decreases very rapidly 
until, in the larger systems, a plateau p + = L~ 2 /2 (i.e., an average 
of 1/2 of active site in the system, due to averaging over up-spin and 
down-spin ground states) is reached. One large collection of positive 
excess is moved around the system in samples with positive magne- 
tization. This packet of excess annihilates broadly distributed sinks 
having a magnitude approximately proportional to A. The density of 
negative-excess sites (sinks) decreases very rapidly, though the rate 
of decrease slows somewhat after the coalescence of positive excess 
into a single packet. 



one sink. For a pure fractal, this would be used to estimate 
the box-counting dimension of the set. (One other charac- 
teristic of the final topography, which we study in Sec.|V] is 
the geometry of paths along unsaturated bonds that start from 
sinks: these are the paths followed in global updates to iden- 
tify reachable sites. At small disorder, these paths have trivial 
geometry.) 
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Figure 5: Plot of the scaled height distribution i 1//2 it 1 ^ 2 P(u) vs. the 
sample-size normalized height, L u, measured at the termination 
of the push-relabel algorithm for one-dimensional samples in the 
limit of weak disorder. The data collapse well to a constant, indicat- 
ing that the distribution of heights behaves as P(u) ~ u^ 1 ^ 2 from 
small scales up to heights u « 10~ 2 L, to within the numerical error 
bars. The error bars represent la estimates for the sample-to-sample 
variation; the errors at a given L are correlated, as the distribution 
at each scaled distance was computed from a single set of samples 
(the number of samples was at least 1700 at each sample size). The 
global update frequency was set at T = 2L. 

1. Remnant sinks for d = 1 

After the terminal global update, the sites i that are not 
sinks have a height label m that equals the lattice distances 
to the nearest sink. Fig.|5]displays a plot of the scaled height 
probability distribution L 1 / 2 it 1 / 2 P(u), where P(u) is defined 
as the probability of a site i to have a height u, for samples 
of various sizes L at small A. The plots show that P(u) is 
well fit by a single power-law behavior, P{u) ~ u~ T , with 
t = 0.500 ± 0.005. (This very small error bar in r is es- 
timated by finding the range of values which give a plateau 
to within statistical error for P(u)L T , for u > 100 and 
2 097152 > L > 16 384 and assuming that the corrections 
to scaling are very small, so that deviations from a plateau 
represent only a statistical error in the exponent.) 

We also characterized the sink distribution using N(b). The 
estimates for this function are plotted in Fig. |6la) for d = 
1. This data is also fit by a single power law, with N(b) ~ 
b-°, where D = 0.500 ± 0.005 (statistical error only). This 
estimate for the fractal dimension is consistent with D = | 
for the set of remnant sinks in the small disorder limit. 

The result for fractal dimension D and height-distribution 
exponent r are related. In one dimension the distribution func- 
tion P(ui) is precisely related to the number n s (£ ss ) of sink- 
to-sink gaps of length i ss via 

oo 

P( Ui ) = n s (2u l -1) + 2V K(2fc) + n s (2k + 1)] . (2) 

k=Ui 

In the continuum height limit, then, taking ^gr^ ~ n s and 



7 



E 1 


1 1 lllllj 1 


1 1 lllllj 1 


1 1 Mill 
















* 


t= 8 














t= 16 












A 


L= 64 












♦ 


L= 256 












□ 


L= 1024 












• 


L = 4096 




Ill 


(b) d=2 










1 


1 I 


1 I 





1 


QMIIII 



12 14 

10 10 10 10" 10 




10° io 1 io 2 io 3 

b 

Figure 6: Plot of normalized box-counting data N(b)L~ d ^ 2 vs. box 
size b at A — » 0. The quantity N(b) counts the number of non- 
overlapping volumes of linear size 6 that intersect the set of remnant 
sinks. The data shown here is averaged over samples for the auxiliary 
field configurations computed at the termination of the standard (non- 
blocking) push-relabel algorithm. 



P ~ J°° n s (£ ss ) d£ S s gives t — D, consistent with our esti- 
mated values. 

The structure of the sinks is very suggestive, when one con- 
siders the apparent relationship between the dynamics of the 
push-relabel algorithm in the limit of small-disorder and stud- 
ies of the A + B — > reaction. In particular, Leyvraz and 
Redner studied the fractal structure of the B particles in the 
limit that the B particles are immobile. The primary differ- 
ence between their analysis and this model is that the A + B 



8 



reaction was considered for diffusive A particles. Here, the A 
particles (corresponding to positive excesses) move directly 
towards the nearest B particle (negative excess sites). 

Our presumption will be that this distinction affects only 
the relationship between length scale and times. In the diffu- 
sive case, the time scale t gives a length scale t ~ t 1 ! 2 . In the 
case of the push-relabel algorithm, we will take t ~ £ln(£) 
or, equivalently up to logarithms of logarithms, £ ~ t(lni) -1 . 
This dependence comes from the linear "attraction" between 
the positive and negative excess sites, with a logarithmic cor- 
rection reflecting the changes in direction that take place upon 
annihilation of positive and negative excess on successively 
larger scales, consistent with the finite-size scaling of the run- 
ning times. 

Given this correspondence, we expect the same domain 
structure of the negative excess sites in the final states as for 
the B sites in the annihilation reaction with immobile B par- 
ticles. Leyvraz and Redner, using random walk arguments 
to sum up the densities of A vs. B particles, find the distri- 
bution of the distances dsB between neighboring B sites to 
have the form P(cIbb) ~ rfgg 2 - Identifying our £ ss with 

d B B gives dN %~^ B ^ ~ n s (£ ss ) ~ £7s , which is numeri- 
cally consistent with our statistics for the final configuration 
from the push-relabel algorithm. It is worthwhile to note 
that this type of picture is also in agreement with the struc- 
ture of the ID RFIM ground state reported by Schroder et al. 
12^1 . which was also derived using absorbing states of ran- 
dom walks. These connections support a unified picture of 
the dynamics of the push-relabel algorithm, the structure of 
the RFIM ground state, and the previously distinct study of 
annihilation reactions. 



2. Structure of remnant sinks, d > 1 

The results we obtain for small disorder in higher dimen- 
sion are somewhat more complex. Scaling behavior is also 
seen, though the results can depend on the details of the algo- 
rithm. In the small disorder limit, the paths to the sinks are still 
linear, i.e., non-fractal and straight (the paths of course must 
be linear at all disorders when d = 1). However, the frac- 
tal structure of the remnant sinks is more apparent in higher 
dimensions. In some cases there appears to be a new length 
scale, intermediate between the microscopic length and the 
sample size, that characterizes the final topography in sam- 
ples with net negative magnetization. 

The topography at the termination of the algorithm is dis- 
played for d = 2 in Fig.Qfor two variants of the PR algorithm. 
One variant is standard: pushes are executed without regard to 
the status of the destination site. We also consider in this sec- 
tion a "blocking" variant of the algorithm that forbids coales- 
cence of positive excess. Whenever a push is attempted in the 
blocking version from a site i with height m to a site j with 
Uj = Ui — 1, the algorithm first checks whether the destination 
site already has positive excess. If there is positive excess at j, 
the algorithm does not push excess in that direction, but also 
does not relabel site i (a relabel would be executed if there 
was no push due to saturated bonds with = 0). This non- 



fa) 



(b) 




Figure 7: [Color online] Topography of the height field upon termi- 
nation of the push-relabel algorithm, for a d — 2 sample of linear 
size L = 512, with weak disorder (A — > limit) for two different 
modifications of the push-relabel algorithm. The sample shown has 
negative magnetization. The lightness of the color (or grayscale) in- 
dicates the magnitude of the height field. Black pixels indicate spins 
with negative excess (height m = 0), dark pixels indicate smaller 
heights, and light regions indicate higher heights. The height field 
and remnant sinks are shown for (a) the push-relabel algorithm with 
no limits on coalescence and (b) a modification where positive excess 
is not allowed to coalesce (the "blocking" version). 



blocking variant slows down the algorithm some, but prevents 
the coalescence of the positive excess into a single packet. The 
number of positive excess sites does not decrease as quickly as 
in the non-blocking version, with the number of excess pack- 
ets remaining comparable to the number of sinks, until N + 
or N~ approach L d / 2 . As can be seen in the figure, the sites 
with negative excess are clearly clustered on several scales for 
both variants of the algorithm. 

We studied this clustering for both variants by again apply- 
ing box counting for the sites with negative excess. The box- 
counting data for the non-blocking variant, with N(b) nor- 



9 



malized by the number of sinks N(0) ~ L d ^ 2 , is plotted as 
L- d l 2 N{b) vs. box size b in Fig. |5Jb,c). We examined this 
data for scaling behavior. The error bars are small enough 
that plotting the local exponent or slope of the ln-ln plot was 
useful in studying the scaling. The scaling ansatz is that there 
is a crossover in the scale-dependent fractal dimension D{b) 
at a scale L x , 



I lllllllj I lllllllj I lllllllj I lllllllj I lllllllj I I III! 



D(b) = D{bL- x ) 1 



(3) 



with D(z) constant at small z and crossing over to D = d 
at large argument. The quantity used to estimate D(b) is 
the discretized logarithmic derivative, 5 \n[N (b)] / 5 \n(b) = 
|ln[iV(2 1 / 2 6)] -ln[iV(2- 1 / 2 6)]}/ln(2), whose negative 
gives an effective dimension when plotted as a function of 
b. The scaled plot of this estimate for D(b) vs. L~ x b is 
shown in Fig. [8] for our best fit values X w 1.0 (d = 1), 
X w 0.55 (d = 2) and X w 0.5 (d = 3). Our estimates 
for the systematic error bars for X are somewhat smaller for 
d = 1 (less than 0.08) than for d = 2,3, where variations 
in X of the order of 0.12 provide plausible, but less clean, 
collapse to a single curve for D. The values of D{b) or 
D(b) at small argument provide an estimate for the effective 
box-counting dimension of the sinks at small scales. The 
d = 1 data, as discussed in the previous section, give a 
fractal dimension consistent with D = I at small scales. The 
d = 2 data are consistent with a convergence to an effective 
dimension of D sa 0.4 ± 0.1 over about 1 decade in scale 
at the smaller scales for the largest samples (4096 2 ) where 
the longest plateau in effective dimension is seen. The d = 3 
data are consistent with D < 0.2 at small scales. At large 
scales, b > L x , the data are quite consistent with D = d, i.e., 
a scale-independent density of sinks. 

The logarithmic derivatives of the box-counting data for the 
blocking variant of FIFO in d — 2, 3 are displayed in Fig.|9] 
(For the case d = 1, we find that while the number of active 
sites evolves differently, with the assumed form 7V + ~ L d l 2 
consistent with the data, the dimensions D w 0.5 and scal- 
ing X sa 1.0 from the box-counting data are the same for the 
blocking and non-blocking variants.) The best collapses are 
seen for X slightly less than 1, X = 0.93, 0.90 for d = 2,3, 
respectively, with an estimated error of about 0.1. At small 
scales in d = 2, the fractal dimension D is approximately 
0.9 ±0.1 at small scales. For d — 3, the convergence is less 
clear, but is consistent with D between 1.0 and 1.5 at small 
scales. Assuming a single fractal dimension over all scales, 
except for corrections near b = L, i.e., taking X = 1, would 
give D = d/2, consistent with our data for the blocking vari- 
ant. 

Comparing the results for the distribution of sinks in these 
two variants of the algorithm and the density of packets dis- 
cussed in Sec. IIVBI leads to a consistent qualitative picture 
for the interacting evolution of the active sites, sinks, and to- 
pography. In the non-blocking variant of the push-relabel al- 
gorithm, it appears that a scale is frozen in once the positive 
excess coalesces. The packet containing the positive excess 
moves all about the sample, cancelling each small negative 
excess site it encounters and changing directions to follow the 
height field Ui, until it is exhausted. This cancellation acts at 
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Figure 8: [Color online] Plot of the local dimension --^^p of the 

set of remnant sinks plotted vs. scaled box size bL~ x . The data 
is sampled at the termination of the standard (non-blocking variant) 
push-relabel algorithm, for d — 1, 2, 3. The sample sizes are L = 
2 8 , 2 10 , . . . , 2 20 from the rightmost to the leftmost curves for d = 
1, L = 2 7 , 2 8 , . . . , 2 12 from the highest to the lowest curves for 
d = 2, and L = 2 3 , 2 4 , . . . , 2 8 for d — 3 (these curves overlap 
significantly). In the color version, the line colors run through the 
sequence dark green, orange, magenta, black, blue, red, light green, 
from smallest to largest L. The dashed light horizontal lines indicate 
local fractal dimension of D — d. 



small scales. The non-blocking variant, however, maintains 
roughly equal numbers of positive excess packets and nega- 
tive excess sinks, so that the positive excess packets have a 
length scale that grows as the length scale of the sinks does, 
throughout the algorithm. The length scale for the topogra- 
phy of the height field is not set by the "freezing-out" of the 
positive excess into a single packet, then, but by the sample 
size. This leads to a scaling for the final distribution of sinks 
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Figure 9: [Color online] Plot of the local dimension — of 

the set of remnant sinks vs. scaled box size bL~ x averaged over 
samples at the termination of the blocking variant of the push-relabel 
algorithm, for d — 2 and d = 3. The sample-size ranges are L = 
2 5 , 2 6 , . . ., 2 12 for d = 2 and L = 2 1 , 2 2 , . . ., 2 8 for d = 3, 
with the leftmost curves corresponding to the largest samples. In the 
color version, the line colors run through the sequence dark green, 
orange, magenta, gray, blue, red, light green, black, from smallest 
to largest L. The error bars are determined by la sample-to-sample 
fluctuations. 



with Iwl, consistent within the large systematic error bars 
for X in d = 2, 3. The non-blocking variant is more consis- 
tent in microscopic rules and with results for the annihilation 
dynamics of particles of constant charge. 



V. GENERAL DISORDER AND CRITICAL SLOWING 
DOWN 

We now consider larger values of A, where the bonds can 
be saturated and block flow, making the rearrangements of 
excess more complicated than the A — > limit. We focus on 
the true transition at A = A c in d = 3, but have also stud- 
ied the topography of itj and timing for A > A^ in d = 1,2. 
Ogielski noted 1 13] that the PR algorithm takes more time to 
find the ground state near the transition in three dimensions 
from the ferromagnetic to paramagnetic phase. For the high- 
est priority queue (HPQ) algorithm, identical to LPQ except 
that active sites with the highest height are kept at the front 



of the list, Middleton and Fisher H 16111911 studied this critical 
slowing down more extensively in the case of d = 3 and work 
has also been carried out for d — 4 [35]. Here, we compare 
the scalings for the timing of FIFO and LPQ algorithms. The 
position of the peak in running time can itself be used to esti- 
mate A c , if one makes the natural assumption that the critical 
slowing down is greatest exactly at A c . This critical slowing 
down is certainly reminiscent of the slowing down seen in lo- 
cal algorithms for statistical mechanics at finite temperature, 
such as Metropolis, and even for cluster algorithms 1 3 Ol 13611. 

Critical slowing down results from the long length scales 
that arise near a continuous phase transition. To connect the 
topography of the w,; to the physical system, we conjecture 
that the scale for the heights m is given by the correlation 
length £. This assumption is natural: the maximal height 
in a domain of linear size £ must have a scale ~ £ for ex- 
cess to be transported across a domain. The spin-spin cor- 
relation functions die off rapidly over longer ranges, so that 
the ground-state computation need not rearrange excess over 
scales greater than £. This conjectured relationship fits the 
numerical data well. 



A. Timing for d= 1,2 

In one- and two-dimensional systems, there is a system- 
size-dependent crossover from large A to small A behavior. 
Above this crossover scale A X (L), physical quantities such as 
correlation length and algorithmic quantities such as Npp_/n 
are expected to converge at large enough L to an infinite- 
volume limit. This limit is expected to exhibit critical diver- 
gences with a critical point at A = 0. 

Our unsealed data for the running time in d = 1 and d = 
2 for FIFO are plotted in Fig. ^3 From the known scaling 
for the correlation length £ ~ A~ 2 1251 l29t 13711 and taking 
Npfi/n oc £ln£, one expects a straight line fit for Npn/n 
when plotted vs. A~ 2 In A. We do not find convergence to a 
single exponent, but instead effective power law ranges from 



(A" 2 lnA) u - 7 over 



N PR /n ~ (A -2 In A) - 25 to N PR /n 

the range A = 0.02 to A = 0.005 for L = 2 24 . It may be that 
even larger samples are needed to see the expected divergence 
in NpR, Despite this, over the range L = 2 14 — > 2 24 , we 
do find that the location of the peak in the running time is 
consistent with the scaling A^ ~ L~ x / 2 . The peak in the 
running time per site therefore does occur when £ ~ L. 

A similar situation exists in d = 2: the divergence of 
the running time with A is not fit well by the simplest scal- 
ing expectations. However, the location of the peak in the 
running time does scale in the expected fashion. The fit- 
ted location of the peak for the running time L~ 2 Np R (A) 
is plotted in Fig. ^] This data is quite consistent with 
L - £ A- 2 y e A o/ A ~ 2 01, with best fit values £ = 17 ± 4, 
A 2 , = 1.3 ± 0.2, and y = 1.1 ± 0.2, citing 2<T-error bars and 
using all of the data points for L > 16. This error estimate 
is consistent with the variation that one gets by changing the 
fit to L > 64. The error bars are large enough that y = 1 is 
certainly an acceptable value. 
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Figure 10: Running times for FIFO algorithm, (a) d — 1 and (b) 
[color online] d — 2 as a function of disorder A and normalized 
by sample volume. The error bars are small for d — 2, so they are 
shown only for the largest size, where the error bars are largest. 



B. Topography for d — 1,2 
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Figure 1 1 : [Color online] A plot of the system size L vs. A^T (cir- 
cles), where A x is the fit to the locations of the peaks in iVpR dis- 
played in Fig. llQI The dashed lines show a fit to the simplest expected 
form L ~ exp f— Ao/A^) for L > 256 while the solid line is a fit 
to all data (L > 16) using the form L ~ A~ 2y exp(-Ao/A 2 ), with 
best fit values £ = 17.3, y = 1.07, and A 2 , = 1.27. 



directly from the fit to the location of the peak in the timing. 
The scaled data shown in Fig. ^] are therefore in agreement 
with a scaling form 

P 2 (u,L,A)~C 1 P2(u/0- 

e A o/ A gave a somewhat worse scaling 



A fit assuming £ 
collapse. 



In contrast with the attempted scaling for the timing 
data, the topographical data for d = 1 at larger A ex- 
hibit a clear scaling with the expected power-law behaviors. 
Fig. presents an example of this scaling for the fraction 
Pl(u, L, A) of sites with height u for d = 1. The plot is of 
the scaled height distribution function A _1 m 1 / 2 Pi(u, L, A) 
vs. the height normalized by the correlation length, uA 2 , for 
various values of L and A. This plot assumes a correlation 
length f (A) ~ A" 2 , with £(A) < L, that P(u) ~ u" 1 / 2 , as 
it is in the limit of small A, and a properly normalized P(u) 
(J du P(u) = 1), which together give a scaling form 



P 1 {u,L 1 A)^Au- 1 ^ 2 P 1 (uA 2 



(4) 



with P(z) constant for small z. The data collapse very well 
for the range 1 > A > A X (L). 

The data for p2(u) also exhibit scaling in d = 2. The 
data then collapse well for disorders A>A X (L), if we 
choose a scale for the maximal u that is proportional to 
A~ 2j, e - ' A °/ A ) oc £, where we take the values of y and Aq 



C. Timing for d = 3 

Our data for the running time iVpR lend themselves well 
to scaling about the critical point. We have not attempted to 
separately infer the critical value A c and the correlation length 
exponent v, which determine the correlation length £ via 



(5) 



but instead use the values determined in Ref. 1 16]. These val- 
ues, A c = 2.270(5) and v = 1.37(4) were found from scaling 
of the stiffness (energy change due to a change in boundary 
conditions) and spin-spin correlation functions and are con- 
sistent with, e.g., the location of the peak in the specific heat. 

There is then one parameter to fit, the dynamic exponent z, 
which describes the divergence in the running time at A c , if 
one assumes the scaling 
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Figure 12: Scaling plot for the distribution of the heights at general 
A for d = 1 and L — 2 14 , 2 18 , and 2 22 , assuming the scaling form 
Eq. J4j. The shape of the symbol indicates the value of A, as indi- 
cated by the shape of the solid symbols in the legend. The color of 
the symbol indicates L, with white-filled symbols for L — 2 14 , gray- 
filled symbols for L = 2 18 , and black-filled symbols for L — 2 22 . 
The data for a broad range of A, including those values that do not 
lie in the scaling range 1 > A > A^ (L) and are therefore not ex- 
pected to collapse, are shown. The points that do not collapse to the 
single curve are indicated by the region labeled A < A x (L) at the 
top of the graph. The points that do not fit for well to the right of the 
collapsed data are for A = 1, as indicated by the legend, where the 
correlation length is comparable to the lattice spacing. The data that 
are in the scaling range give a good collapse to a single curve, with 
many points overlapping to the extent that they are not visible. The 
scaling curve is flat over the range 10" 4 < A-V: < 10 _1 . 



Figure 13: Scaling plot for the distribution P(ui) of the final heights 
Ui at general random-field strength A for d — 2, with L == 2048 
and L = 4096. The correlation length £ which sets the scale for the 
heights Ui is taken to diverge as £ ~ A~ 2y exp(Ag/A 2 ), with the 
parameters y = 1.07 and Ao = 1.27 taken from the scaling for the 
peak running time CFig. ll li . Curves and points that are expected to 
be outside of the scaling range are included for comparison. Besides 
small values of m (L — 2048 points near the top-center of the plot), 
the data that does not collapse is where A > 1.25, where £ is becom- 
ing comparable to the lattice size, and for L = 2048, A = 0.55, 
which is on the low- A side of the running time peak (see Fig. 1 101 . 
The curves over the range 0.55 < A < 1.00 collapse relatively 
well; by the locations of the finite-size peaks in the running times of 
Fig.Qol this range of A corresponds to almost two decades in length 
scale: 64 < L < 4096. A fit to the simpler £ ~ exp(Ao/A 2 ) with 
a fit parameter Ao gives a worse scaling collapse. 



for the number of PR steps per site, where w(x) ~ x~ z at 
large x and w(x) ~ \x\~ z ln( ) as x — > — oc, to be con- 
sistent with convergence to constant L~ d NpR at A > A c 
and L~~ d NpR ~ ln(L) for small A (note that the coefficient 
of this logarithm in the limit of large L is probably different 
from what we find here, given that the true logarithmic be- 
havior may not have been reached, as discussed in Sec. llV Al l. 
Our best fits to this scaling form are plotted in Fig.^Jfor the 
LPQ and FIFO data structures. Our estimates for the dynamic 
critical exponent z are distinct for these two structures, with 

zlpq = 0.93 ± 0.06; z F i FO = 0.43 ± 0.06 . (7) 

The error estimates reflect the range of fits that are consistent 
with a correction to scaling that is not too large; the statisti- 
cal error bars are quite small for this data. The convergence 
for the FIFO data may be slower as the data diverge with L 
more slowly, so that corrections to scaling arising from, for 
example, constants are more evident. 

Given the scaling in the data about the thermodynamic criti- 
cal point, it is natural to attempt to explain the critical slowing 
down as being due to the divergence of a correlation length 
£. The finite-size effect (scaling with L) then reflects how the 
running time diverges with £ in the infinite volume limit. The 
difference in the scaling of the running times for the two dif- 
ferent data structures indicates different scaling with respect 



to £. This difference reflects the order in which the operations 
are carried out. For the LPQ algorithm, only a subset of the 
sites, those with the lowest Ui, are subject to PR operations. 
This is in contrast with FIFO, where all active sites are con- 
sidered cyclically. The FIFO structure leads to a more rapid 
coalescence of active sites and cancellation of positive packets 
with sinks. It appears that it then takes fewer PR operations 
to transport the same quantity of excess across a domain as 
the domain size increases. In the LPQ implementation, the 
value z « 1.0 suggests that order £ operations are carried out 
on average in domains of si ze j d (note that z « 1.0 in 3D 
for highest priority queues V ila. B5ll . This is consistent with 
each site being relabeled an average of a multiple of £ times. 
The LPQ algorithm coalesces the positive excess by sweep- 
ing across the height. The FIFO algorithm must relabel the 
sites in such a way that only a small fraction of the sites are 
relabeled to height Ui ~ £, with P(u) ~ u z ~~ 2 if P(u) is de- 
scribed by a simple power law and the number of relabels is 
proportional to iVpR. Nonetheless, taking £ as a cutoff for the 
distribution of Ui in FIFO is consistent with the assumption 
that some portion of the excess must be rearranged across the 
width of a domain. 

We note that there are many distinct scaling plots that can 
be made from the timing data. The running times for posi- 
tively magnetized and negatively magnetized samples differs 
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significantly for small A, due to the asymmetry between the 
negative excess sinks and the positive excess packets at ac- 
tive sites. For large A, the algorithm regains its symmetry, 
in that the average magnetization is near zero and any sam- 
ple is a nearly even mix of up-spin and down-spin domains. 
One of the plots that we found to be a sensitive test of scal- 
ing was to compare the running time A/pR between samples 
with opposite magnetizations. In particular, the relative rms 
fluctuations <7^ pR and cr^ pn in the operation count for sam- 
ples with positive and negative magnetization, respectively, 
vary in a rapid fashion near criticality. The dimensionless ra- 
tio = crtr I cr7 T is a non-monotonic function charac- 

iVpR / iVpR 

terizing the distribution of running times. This function of A 
and L collapses very well for L > 16 with no adjustable pa- 
rameters (taking A c and v fixed as stated earlier). A plot of 
this collapse is shown in Fig.^lf or the FIFO algorithm. The 
parameter-free fit (subject to using published values for A c 
and v) provides a clear confirmation of scaling for the distri- 
bution of running times. 



D. Topography at A = A c for d = 3 

The evolution of the auxiliary fields are closely connected 
with the timing of the algorithm, as seen in Sec. ll Vl for A — > 0. 
We have studied the final topography of the height field and 
the paths that lead to the sinks for d — 3 and general A. We 
use this information to study the characteristic topography at 
A near A c when using the FIFO data structure. 

A study of the box-counting data for the sinks gives less 
information than for A = 0. For any A < A c , there will 
be a finite density of minority spins (spins opposite to the net 
magnetization). This density drops very rapidly in d — 3 as A 
decreases below A c , but these minority spins result in a fractal 
dimension of D = 3 at all scales £ greater than the typical 
separation between the minority spins. We found no obvious 
singularity in N(b) for A « A c : the logarithmic derivative for 
N (b) converges rapidly to rj 3.0 for b > 8, for the disorder 
range 2 < A < 3 . The spatial distribution of the sinks is non- 
fractal for A « A c (at least from simplest the box-counting 
perspective; there will of course be some singularities in the 
density and correlations of down spins at A c ). 

A feature of the topography that is much more sensitive to 
the phase transition is the distribution of the u. The height 
values Ui reflect long range correlations of the spins (minor- 
ity sinks give D = 3 but are localized in their effect on the 
Ui). In Fig. ^] we compare the distributions P(u) for the two 
largest systems we studied at finite A, L = 128 and L — 256 
to indicate how L and A appear to affect the height distri- 
butions. At large A, the height distribution decreases expo- 
nentially with A, P(u) ~ e~ u /^ A \ This is consistent with 
spin-spin correlations decreasing exponentially over a charac- 
teristic length scale £(A) in the paramagnetic phase. For A 
significantly above A c , the P(u) are relatively independent 
of L. For A <C A c , the height distribution peaks at a scale 
■u^ ~ L x that grows with system size. This is consistent with 
the data of Fig. [8] which gives a crossover in the spatial distri- 
bution of the sinks (for the non-blocking variant). Above this 



4 
3.5 

3 

,2.5 
2 
1.5 
1 

0.5 




o 


L= 


8 


♦ 


L= 


= 16 


□ 


L= 


32 


■ 


L= 


64 


o 


L= 


128 


• 


L= 


256 




FIFO, d=3,z = 0.43 



2 
(A-AJL 



l/v 



1 

0.8 
0.6 
0.4 
0.2 




o 


L=8 


♦ 


L=16 


□ 


L=32 


■ 


L=64 


o 


L=128 



□ 



(b) 



of 



if 



o 



o 



. o 
♦ o 

o oi 



LPQ, d = 3, z = 0.93 



2 4 6 
(A-AU" 1/V 



8 10 12 



Figure 14: [Color online] Plot of scaled running time in d — 3 for 
(a) the FIFO data structure, with a fit value of z = 0.43, and (b) the 
LPQ data structure, with dynamic exponent z — 0.93. The scaling 
fit assumes A c = 2.270 and v — 1.37 as fixed parameters. 



length scale ~ L 5 , the distribution of sinks becomes uni- 
form, leading to a cutoff in the itj, or distance to the nearest 
sink, above that scale. Near criticality, the distribution P{u) 
must crossover between these two behaviors, one decreasing 
at small u and the other increasing. The distribution changes 
rapidly with A and L in this region. We find that the criti- 
cal distribution is not varying exponentially rapidly for some 
A«A C . 

One other feature of the topography that we have investi- 
gated is the structure of the paths that connect down-spin sites 
to a sink at the final step of the algorithm. When A is small, 
these paths must be linear, as there are no saturated bonds. 
Near criticality, many bonds are saturated and there is the po- 
tential that the paths to a sink could be rather torturous, as they 
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Figure 15: [Color online] Plot of the dimensionless ratio R+- of 
the rms fluctuations o"^ pR in iVpR, (time to find the ground state) 
positive magnetization samples to the corresponding quantity o"jv pR 
in negative magnetization samples plotted vs. the scaled disorder 
(A - A C )L 1// ", for samples of sizes L = 8, 16, 32, 64, 128, and 
256. The values A c = 2.270 and v = 1.37 were assumed, based 
on scaling of physical (not algorithmic) quantities, so this parameter- 
free fit is a direct check of scaling for the distributions of running 
times. 



would need to avoid saturated links. Numerical studies were 
carried out to measure the paths to the sinks. When carrying 
out a global update, we maintain a destination label that gives 
the location of the sink to which any positive excess will flow. 
For fixed height of each site with u < oo, the Euclidean dis- 
tance r of that site to the sink given by the destination label 
is computed. The FIFO data can be used to estimate a fractal 
dimension df,u ~ r d f , for the paths which guide the pushing 
of the excess. Near A c , it appears that the paths are nearly 
linear. For u > 10 and A < 2.30, the Euclidean distance r 
was linear in height u, with an effective fractal exponent in the 
range 0.9 < df < 1.1 for samples of size L = 256. 

VI. SUMMARY 

In this paper, we have studied the dynamics of the auxil- 
iary fields used in the push-relabel algorithm to find the exact 
ground state of the random-field Ising magnet, a prototypical 
glassy model. These dynamics, and the final state of the fields 
found during the solution, reflect both the underlying physi- 
cal model and the special dynamics of the algorithm. There 
is some freedom of choice for these dynamics: we studied 
primarily the FIFO queue with global updates, which is the 
fastest algorithm we have used and is also most directly sim- 
ilar to the synchronous evolution used to model physical dy- 
namics. The evolution of the auxiliary fields is roughly a re- 
arrangement and cancellation of the locally varying external 
fields. The extent of this cancellation determines the size of 
the domains in the RFIM ground state. 

In the limit of small disorder, the dynamics of the field 
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Figure 16: Plot of the height probability distribution P(u) for differ- 
ent values of disorder A and size L, for d = 3 near the critical value 
A c = 2.27. The symbols and error bars are for a sampling of points; 
the solid lines are measured for every integer value of u. (a) At low 
A < Ac, the distribution of heights it has a peak at an Independent 
scale, (b) Near A c , the distribution of u is much less dependent on 
u and is nearly constant over a large range of u (on this scale, power 
laws are nearly constant), though the form of P(u) varies rapidly for 
A ~ A c . (c)High Avalues, here A = 2.40 > A c give a distribution 
P(u) ~ e"" /e . 



of excesses is that of a potential-driven annihilation reaction, 
A + B — > 0. The algorithmic time per spin to find the ground 
state for N spins is apparently linear in ln(iV), though a sim- 
ple logarithmic behavior is seen only at large system sizes. 
The densities of the mobile positive excess fields (active sites) 
and the immobile negative excess fields (sinks) decreases very 
rapidly with time. When the positive excess is allowed to co- 
alesce, the number of positive sites becomes of order unity 
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before the ground state is found. This coalescence leads to 
a cutoff in length scale L x in a sample of linear size L that 
limits the range of scales over which the sinks at the end of 
the algorithm can be described as fractal. If coalescence is 
forbidden, the remnant sinks appear be described by a fractal 
measure all the way up to scale that is nearly equal to L. The 
fractal dimension in d = 1, D » 0.50, is consistent with ran- 
dom walk arguments that have been developed for diffusive 
annihilation-reaction dynamics l2lll . For d = 2, D w 0.4 and 
the fractal dimension D < 0.2 for d = 3. 

Diverging length scales in the RFIM influence the running 
time for the algorithm at general disorder values. While the 
running time itself is not easily scaled, the peak location for 
d = 1,2, A X (L), scales as expected from calculations for the 
divergence of the correlation length £ with A in the RFIM 
J23I I24I l39ll . In particular, for d — 2, the location of the 
peak in the running time for a finite system and the distribu- 
tion of the height fields near the transition give evidence for 
the predicted form for the power law corrections to the expo- 
nential dependence of £ on A -2 . Further examination of the 
data may be useful in explaining how Np~r(A) diverges. In 
d — 3, the scaling of the running time (and its distributions, 
as seen in Fig.1151 are quite consistent with values [16] for the 
critical value A c and correlation length exponent v obtained 
from simulations for the physical properties of the RFIM. For 
A w A c in d = 3, the distribution of the field corresponding 
to the potential Ui generated by the sinks is nearly constant in 
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Ui. For small A > A x , P(ui) <x u i for u <C £ in d = 1 
and varies very slowly for itj up to £ in d = 2. 

We believe that these consistent and strong connections be- 
tween algorithm dynamics, the physics of the RFIM, and the 
mathematics used to describe annihilation processes provide 
insight into both the push-relabel algorithm and the RFIM. 
One is tempted to speculate that relating the rearrangement 
of fields in the push-relabel algorithm to finding the domains 
in the RFIM might help in analytic approaches, either in fi- 
nite dimensions or on hierarchical lattices |40]. The running 
time of the PR algorithm is directly related to the average of 
the potential or height field; distinct algorithms give different 
scaling for this average, though the physical ground state that 
is found is identical. A better understanding of this dynamical 
exponent z given by the average height, is likely to shed light 
on both the algorithm and the RFIM ground state at criticality. 
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